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ABSTRACT 

Mergers of gas-rich galaxies are key events in the hierarchical built-up of cosmic 
structures, and can lead to the formation of massive black hole binaries. By means 
of high-resolution hydro dynamical simulations we consider the late stages of a gas- 
rich major merger, detailing the dynamics of two circumnuclear discs, and of the 
hosted massive black holes during their pairing phase. During the merger gas clumps 
with masses of a fraction of the black hole mass form because of fragmentation. Such 
high-density gas is very effective in forming stars, and the most massive clumps can 
substantially perturb the black hole orbits. After ^10 Myr from the start of the merger 
a gravitationally bound black hole binary forms at a separation of a few parsecs, and 
soon after, the separation falls below our resolution limit of 0.39 pc. At the time of 
binary formation the original discs are almost completely disrupted because of SNa 
feedback, while on pc scales the residual gas settles in a circumbinary disc with mass 
^ 10^ Mq. We also test that binary dynamics is robust against the details of the SNa 
feedback employed in the simulations, while gas dynamics is not. We finally highlight 
the importance of the SNa time-scale on our results. 

Key words: black hole physics - hydrodynamics - galaxies: formation - galaxies: 
evolution - galaxies: nuclei 


1 INTRODUCTION 


Understanding the black hole binary formation path in 
galaxy’s mergers is a key step towards understanding the 
mode of assembly of massive black holes (MBHs) across cos¬ 
mic history ( Di Matteo, Springel Hernquist||2Q05 Sesana 


et al.|2Q14 Dubois, Volonteri Silk|2014[ ). The masses and 

spins of MBHs are sculpted by complex processes that in¬ 
volve episodes of accretion and mergers, during the hierar¬ 
chical growth of structures. The coalescence in a single, more 
massive BH (induced by gravitational wave emission) is ex¬ 
pected to occur whenever the two MBHs present in each of 
the interacting galaxies reach sub-parsec separation under 
the action of stellar and gas torques (Amaro-Seoane et al. 


2013). 


The dynamics of MBHs is complex as it occurs in the 
time-varying environment of a galaxy merger where stars, 
and cool gas, turning into new stars, evolve on comparable 
time-scales (see |Colpi 2014| for a review), and is custom¬ 
arily described as a three-step process: (I) pairing during 
the merger of the galaxies embedded in their dark matter 
haloes, resulting in the formation of a Keplerian binary, (H) 


migration by gas or/and hardening by stars , and lastly (HI) 
gravitational wave driven inspiral. 

The seminal paper by [Begelman, Blandford Rees| 
(1980) showed that in an already established spherical 
galaxy remnant devoid of gas, the fastest of these three 
phases is that of pairing under the action of dynamical fric¬ 
tion against stars. The MBHs are driven down to parsec sep¬ 
arations and form a Keplerian binary, when the stellar mass 
enclosed in their orbit drops below their mass. IBegelman^ 


Blandford & Rees (1980) further noticed that hardening of 


the binary by encounters with individual stars in phase (H) 
is the most critical and long-lived phase, representing a bot¬ 
tleneck to the path to coalescence, at least for MBHs with 


masses larger than ~ 2x lO^M© (see e.g. Merritt & Milosavl- 


jevic 


2005 Merritt, Mikkola & Szell 2007). The so-called 


last parsec problem for the most massive MBHs has been 
recently allevi ated after considering more realistic galaxy’s 
background^ [Begelman, Blandford Rees (1980) mention 


^ We will not consider here MBH pairing in pure stellar environ- 
ments and refer to|Preto et al. ](|2011| ); |Khan et ^(|2Q13| ) ;|Vasiliev,| 
[AntoninT Sz Merritt] ( |2014[ |2015| ) "and [Sesana Sz Khan| ( |2015| ) for 
recent findings. 
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on the possibility that gas can accelerate the orbital decay. 
In this case, torques from gas distributed around the binary 
in the form of a circumbinary disc are central for driving the 
binary down to coalescence but they depend on how mas¬ 
sive and long-lived are these discs and little is known about 
their formation, structure and lifetime in the relic galaxy 


(Cuadra et al. 2009 Roedig et al. 2011 2012 del Valle & 
Escala||2012 ). 

Thanks to major advances in numerical computing, a 
number of studies have begun tracing phase (I) of MBH 
pairing in the very early stages of a merger when the two 
dark matter haloes first touch, and later when the galactic 
discs (stellar and gaseous) and bulges interact to drive and 
terminate the morphological transformation of the galaxies 
(Mayer et al. 2007| [Roskar et al. 2015). A leap in under¬ 
standing the role of gas during the pairing phase was taken 
when studying minor mergers, i.e. mergers with nominal 1:4 
mass ratios and less. It was demonstrated that thanks to its 
dissipative action, gas deepens the gravitational potential of 
the less massive galaxy reducing the action of tides by the 


primary (Kazantzidis et al. 2005). This prevents the wan¬ 
dering to the lighter MBH in the periphery of the primary 
when sufficient gas is present in the host galaxy, raising the 
question as to whether minor mergers lead in general to 
binary formation (Callegari et al. 2009 2011). Follow-up 


studies have indicated that the dividing line from success 
and failure in forming an MBH binary (MBHB hereon) is 
around 1:10 mass ratios (Bellovary et al. pMol 


Van Wassen- 


hove et al. 2012), but still depends on details such as the 


encounter geometry and gas content. 

Major mergers among gas-rich galaxies represent a nat¬ 
ural path for MBH pairing and binary formation. The or¬ 
bital braking is in these cases driven by gas-dynamical fric¬ 
tion which is faster than dynamical friction fr om stars (|Es- 
cala et ar]|2005| |Dotti, Colpi &: Haardt||200^ [Mayer et al. 


2007||Chapon, Mayer V Teyssier||2013 ). Massive inflows of 

gas during the merger lead to the formation of a circumnu- 
clear disc comprising most of the gas initially hosted in the 
individual discs. It is in this new disc that the two MBH 
sink forming a Keplerian system down to pc scale on a time 
typically of ~ 5 Myr after completion of the merger. A crit¬ 
ical moment was proven to be the time when the two discs 
start colliding and developing shocks due to orbit crossings 
which cancel/redistribute the gas angular momentum. This 
process triggers inflows that lead to the formation of a tur¬ 
bulent Toomre stable, massive nuclear disc (Chapon, Mayer 


& Teyssier 2013). The limit and drawback of these high- 


resolution simulations is that gas is treated as a single phase 
medium described by a polytropic equation of state (with 
index 7/5) which mimics the thermodynamics of a generic 
star-forming region. 


Capelo et al. (2015) studied the large-scale dynamics of 


MBH in a variety of mergers with mass ratio 1:1 down to 1:10 
to explore the black hole accretion history and their dynam¬ 
ics during the pairing phase, down to separations of several 
parsecs. Their focus was mainly in exploring the possibil¬ 
ity of triggering ‘dual’ AGN activity along the course of the 
merger. Only recently state-of-the-art simulations of major 
mergers have achieved enough resolution to detail the MBH 
dynamics down parsec scales, in presence of a multiphase 
gas ( Roskar et al.||2015 ). 

Three lines of studies have been developed in the past. 


showing that the MBH dynamics in star-forming regions is 
more complex than predicted by single-phase models and 
that the black hole environment is far richer than expected. 
The three approaches consist of: (a) simulations of major 
mergers of gas-rich disc galaxies starting from large scales 
(> 100 kpc; Roskar et al. 2015); (b) simulations of two 


gaseous discs of mass ~ 10^ M© colliding on scales of ~ 500 
pc, representing a zoom-in version of very gas-rich galaxy’s 


discs in the verge of merging ( Lupi, Haardt V Dotti||2015 ); 


(c) simulations of a single gaseous disc unstable to frag¬ 
mentation as it forms as a consequence of major gas-rich 
mergers. This last kind of numerical experiment focused on 


2013 del Valle et alj 

|2015) as well as on lighter (< 10 ^M( 

s) 

circumbinary discs (/ 

^maro-Seoane, Brem V Guadra|2013 

2 


With the smoothed particle hydrodynamics (SPH) code 
GASOLINE, case (a) was explored by Roskar et al. ]( |2015| who 
simulated a prograde in-plane 1:1 merger of two late-type 
galaxies (Milky Way like) each having a 10% fraction of 
gas in their discs. Star formation and feedback were turned 
on at the onset of the simulation to generate a multiphase 
medium. 5 Gyr after the start of the simulation, particle¬ 
splitting of the baryonic particles was performed in an ex¬ 
cised zoom-in region to follow the last 100 Myr of evolution 
when the galaxy’s cores touch on scales of ~ 5 kpc. The 
starburst triggered by the merger and stellar feedback cre¬ 
ate a nuclear region temporarily devoid of gas and only after 
~ 10 Myr a nuclear disc is rebuilt. The nature of the multi¬ 
phase gas which develops clumps affects the MBH dynamics. 
The MBHs undergo gravitational encounters with massive 
gas clouds and stochastic torquing by both clouds and spiral 
modes in the disc relents the pairing process. The MBHs are 
kicked out of the plane due to their interaction with clumps 
and this delays the time of binary formation which now is 
~ 100 Myr, about two orders of magnitude longer than in 
the idealised mergers with one-component gas. This finding 


is in line with that found by Fiacconi et al. (2013) who sim¬ 


ulated the dynamics of MBHs in a massive gaseous disc of 
< 10®“® M© subject to rapid cooling and fragmentation, in 
a suite of runs of type (c). The stochastic behaviour of the 
MBH orbit was found to emerge when the ratio between the 
clump mass and black hole mass is greater than unity, as in 
this case the impulsive interaction can modify the MBH or¬ 
bit substantially. Decay time-scales are found to range from 
1 to ~ 50 Myr, under the conditions explored. 

Another type (c) investigation has been performed by 
del Valle et al. (2015), who recently simulated with an SPH 


code the sinking of two MBHs in a massive 10® M© circum- 
nuclear disc with gas forming stars. The orbits of the MBHs 
are erratically perturbed by the gravitational interaction 
with the clumps that form as a result of disc’s fragmen¬ 
tation, delaying the orbital decay of the MBHs if compared 
with similar runs with a one-component gas: typical decay 
times are found close to 10 Myr, when the MBHs are seeded 
in the disc initially at ~ 200 — 100 pc scales. The key result 
which emerges from these new studies is that the MBH dy¬ 
namics is sensitive not only to the time varying gravitational 


^ Noticeably, the pairing of the MBHB in this last study has been 
followed using a high accuracy direct A-body code to evolve the 
distribution of stars previously obtained from an SPH run. 
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background of a merger itself, but also on how fragmentation 
of gas clouds, star formation and supernova (SNa) feedback 
shape and change the thermodynamical state of the gas, 
considered to play a key role in guiding the orbital decay of 
the MBH. The mass distribution of the star-forming clumps 
appears to be a relevant parameter which affects the degree 
of stochastic forcing of the MBH orbit and the distribution 
of the sinking times from ~ 100 pc scale down to 0.1 pc. 

In this paper we simulate the evolution of two gaseous 
discs [type (b)] and of their embedded black holes in the 
aim at studying the MBH dynamics within a multiphase 
gas shaped by cooling, star formation and SNa feedback, 
using the adaptive mesh refinement code (AMR) RAMSES 
( Teyssier]|2002 ) modified to accurately track the MBH dy¬ 
namics as presented in Lupi, Haardt &: Dotti| ( |2015| ) . AMR 
codes like RAMSES are known to better resolve gas shocks 
with respect to SPH codes ( Agertz et al.||2007 ), thus allow¬ 
ing a more accurate description of the gas dynamics when 
the two gaseous discs collide. A key question to pose is: in 
a star-forming medium what is the role of SNa feedback in 
shaping the gas mass distribution around the MBHs? The 
transit from the binary phase H to phase HI of gravitational 
wave inspiral depends on the strength of gas-driven migra¬ 
tion in a circumbinary disc surrounding the MBH binary 


( Cuadra et al.|200^ Shi et al.|2012 Roedig et al.|201lj 2012 


del Valle Escala||2012 ). Here we first attempt to explore 
under which conditions a circumbinary disc forms around 
the two MBHs and how this depends on the recipes adopted 
to model the physics of star-forming regions. 

The paper is organised as follows. In Section 2 we de¬ 
scribe the numerical setup used to simulate the massive gas 
discs embedded in a stellar bulge and the recipes introduced 
to model star formation and SNa feedback. In Section 3 we 
present the results focusing on both the MBH dynamics and 
the properties of the multiphase gas surrounding the MBH 
binary. Section 4 contains our conclusions. 


2 NUMERICAL SETUP AND INITIAL 
CONDITIONS 


We consider the late stages of a galaxy gas rich merger, 
in which both galaxies host an MBH surrounded by a cir- 
cumnuclear disc. Full details of the initial conditions can 
be found in Paper I. Here we simply summarise the basic 
features. 

We initially set each of the two merging nuclei in dy¬ 
namical equilibrium, assuming they are constituted by three 
different components: 


• a stellar spherical structure (termed ‘nucleus’ hereafter) 
described by an Hernquist profile ( Hernquist|I99Q| ), defined 
in spherical coordinates as 


Pb{r) 


Mb a 
27 r r (r + a)^ ’ 


( 1 ) 


where pb{r) is the density as a function of radius r, Mb = 
2 X lO^M© the total nucleus mass, and a = 100 pc the 
nucleus scale radius containing 1/4 of Mb; 

• an exponential gaseous disc with surface density profile 
defined (in cylindrical coordinates) as 


^d(R) 


Md 

2nRl 


exp{-R/Rd), 


( 2 ) 


where R is the disc radius, Rd = 50 pc the disc scale radius 
containing 0.26 of the total disc mass Md = IO^Mq; 

• an MBH with mass Mbh = IO^Mq , at rest in the centre 
of the disc. 


We build the two equal mass corotating gaseous discs, 
each described by 10^ particles, by means of the publicly 
available code gd_basic (see |Lupi, Haardt Dotti|20I5 for 
the algorithm description) and we relaxed them for about 
10 Myr to ensure stability. The discs are initially set at 300 
pc on an elliptical orbit with eccentricity e = 0.3, and with 
orbital angular momentum antiparallel to the angular mo¬ 
mentum of the discs. We stress that each galaxy disc plane is 
in principle uncorrelated to the orbital plane of the merger, 
and, to the first order, the same is valid for the CND^ We 
arbitrarily chose the geometry that maximises the impact of 
the two discs along their orbit and that ensures the highest 
cancellation of angular momentum, enhancing the inflows 
towards the centremost regions. Such a geometry has not 
been explored in the literature yet. The initial conditions 
for the AMR runs are obtained by mapping the gas parti¬ 
cle distribution on the grid using the publicly available code 

tipgridQ 

We perform a total of six simulations using the AMR 
code RAMSES ( Teyssier|2002 ), where we change the prescrip¬ 
tions regarding gas cooling, stellar mass formation (SF) and 
SNa feedback to survey how the implementation of sub¬ 
grid physics affects the evolution of the system. The gas has 
primordial composition, it is optically thin and cools down 
under lines and continuum emission. The maximum spatial 
resolution (at the highest refinement level) for all our sim¬ 
ulations is ~ 0.39 pc and the mass resolution for particles 
forming the stabilising stellar nucleus is 2 x lO^M©. The 
Jeans length is always resolved with at least 4 cells (14 in 
the highest refinement level) and we also add a pressure sup¬ 
port term, modelled as a polytrope with 7 = 5/3 and tem¬ 
perature 2 X 10^ K at the star formation threshold, in order 
to avoid the formation of spurious clumps due to resolution 
limits. 

In order to achieve the best possible treatment of MBH 
dynamics, we adopt the additional refinement criterion de¬ 
scribed in Lupi, Haardt & Dotti ( 2QI5| ). We allow stellar 
particle creation when gas matches two criteria: (i) the gas 
temperature drops below 2 x 10^ K, and ii), the gas density 
in a cell exceeds a pre-defined value. We assume two different 
fiducial values for the SF density threshold, i.e., nu — 2x 10^ 
and nn = 2 X 10® cm“^, where nn is the local hydrogen num¬ 
ber density, and a typical SF time-scale of I.O Myr. The re¬ 
sulting average mass of stellar particles is of ~ 300 M©. Such 
value is significantly more massive than, e.g., what employed 


Amaro-Seoane, Brem & Cuadra (2013), who however sim¬ 


ulated a lighter and more compact system. We checked that 
our prescription results in a gas-to-stellar mass conversion 
rate not lower than the local Kennicutt-Schmidt law. 


^ Here we are neglecting the possible tidal effect exerted by one 
disc on to the other. This effect would tend to align (or antialign) 
the two discs, enhancing the chances of having an orientation 
between the two CNDs similar to the one we assumed as initial 
conditions. 

^ The code is available at http://www.astrosim.net/code/doku. 
php?id=home:code:analysistools:misctools 
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Run 

nil 

(cm“^) 

(Myr) 

Feedback 

ThFBl 

2 X 10® 

10.0 

Thermal 

ThFBh 

2 X 10® 

10.0 

Thermal 

BWFBl 

2 X 10® 

10.0 

Blast wave 

BWFBh 

2 X 10® 

10.0 

Blast wave 

ThFBh.prompt 

2 X 10® 

0.0 

Thermal 

BWFBh.prompt 

2 X 10® 

0.0 

Blast wave 


Table 1. The complete suite of runs. The second column shows 
the density threshold for SF, the third column the lifetime of mas¬ 
sive stars and the fourth column the type of feedback employed. 


In order to model SNa explosions, we consider each stel¬ 
lar particle as a stellar population following a Salpeter IMF, 
and a SNa yield of 15%. We further employ two different 
recipes for the SNa thermal feedback. In both SNa feedback 
recipes the energy budget associated (lO^^erg/M©) is com¬ 
pletely released in the parent cell as purely thermal energy. 
The first criterion (termed ‘ThFB’ after thermal feedback) 
assumes that the heated gas starts cooling right after the 
SNa event; the second criterion (termed ‘BWFB’ after blast 
wave feedback) assumes instead that the energy released by 
SNe is decoupled from the gas radiative cooling, i.e., it is not 
radiated away for ~ 20 Myr ( Teyssier et al.|[20T3 ) and this 
triggers the formation of a momentum-driven blast wave. 
This latter scheme is aimed at modelling non-thermal pro¬ 
cesses energising the blast wave, which are characterised by 


time-scales longer than thermal processes (see e.g. EnBlin 
et al.|[2007 ). We usually assume that no star formation oc¬ 


curs within the two discs before the merger, and that SNe 
explode after a time AtsN = 10 Myr. Stellar mass parti¬ 
cles forming the stabilising bulge are not allowed to release 
energy as SNe. 

We run two further simulations at the highest density 
threshold for star formation (termed ‘ThFBh_prompt’ and 
‘BWFBh_prompt’) assuming no time-lag between star for¬ 
mation and SN explosion. The aim of these runs is to test 
the effects on the global (gas and BHs) dynamics of a max¬ 
imally fast SN feedback, comparing the results to the stan¬ 
dard Ats = 10 Myr case. While simulations with standard 
delay are meant to model star formation as triggered by 
the merger of the two circumnuclear discs, the 0-lag case 
may represent a situation where sustained star formation is 
already in progress at the time of the merger. 

Finally, in order to avoid inaccurate integration of the 
orbits, all runs are stopped when the MBH separation is 
approximatively three to four times the cell length. Table 
summarises our six simulations with the parameter used. 


3 RESULTS 

3.1 Black hole dynamics 

We start by describing the MBH dynamics for the two simu¬ 
lations characterised by standard thermal SNa feedback and 
a typical time for SNa explosions of 10 Myr (runs ‘ThFBF 
and ‘ThFBh’). These two runs are meant to represent a case 
where star formation is indeed triggered by the merger event, 
while gas thermodynamics is governed by standard thermal 


processes. The two different density thresholds for SF are 
used to assess the effects that the efficiency of gas conver¬ 
sion into stars has on the MBH dynamical evolution. 

In Fig. [^we show the MBH projected orbits (left-hand 
panel), and MBH separation versus time (right-hand panel). 
The two MBHs exhibit a peculiar orbital motion, which can 
be explained when considering the gravitational interactions 
between the MBHs and massive gas/star clumps forming 
in the merging discs. Such interactions typically accelerate 
the orbital decay of the MBHs, and a gravitationally bound 
MBHB forms after ~ 10 Myr (the binary formation time is 
indicated as a blue dot in the right-hand panel). In the case 
of ThFBh run (dashed red lines), the MBH orbits appear 
more perturbed, and the orbital decay is somewhat faster. 

Fragmentation of gas, occurring just after the simula¬ 
tion starts, tends to form massive gas clumps, especially 
in the high-density regions surrounding the two MBHs. In 
high-density clumps star formation is very effective, and 
overall, a large fraction of the initial disc gas is converted 
into stellar mass within 10 Myr. This is apparent from Fig. 
[^(left-hand panel), where the stellar mass and the residual 
gas mass are shown as a function of time. The right-hand 
panel of Fig. instead, shows the star formation rate ver¬ 
sus time. A fast increase of star formation occurs initially 
since gas shocked during the disc collision fragments into 
small clumps which immediately convert into stellar parti¬ 
cles. After ~ 2 Myr, only low-density gas survives. Hence, 
star formation is no longer efficient and almost steadily de¬ 
creases in time. In Fig. we plot the mass-weighted gas 
density map at time t = 2.1 Myr, defined as the time of the 
peak in star formation rate (see Fig.[^ right-hand panel). 

In order to quantify the impact of gas clumps on MBH 
dynamics, we estimate the total mass in gas/star clumps, 
along with the clump mass distribution. We define ‘clumps’ 
those gravitationally bound regions that feature a single 
peak in the 3D density field. Fig. shows the total mass 
in clumps as a function of time for run ‘ThFBF. Two dis¬ 
tinct phases can be observed, with a peak in the total mass 
of clumps occurring after a time tpeak ~ 2.5 Myr. The ini¬ 
tial fast growth of the gas locked in clumps is the result of 
the collision between the two unperturbed gaseous discs. In¬ 
deed, gas fragmentation is promoted along the shock surface 
(resulting also in the peak of star formation rate, see Fig.[^ 
right-hand panel). 

Fig. E shows, for the same ‘ThFBF run, the mass dis¬ 
tribution of clumps at four different selected times marked 
as red dots in Fig. We selected two times correspond¬ 
ing to a relatively low total clump mass (~ 1.8 x 10^ M©, 
left-hand panel) and two times corresponding to a larger 
mass value (~ 5 x 10^ M©, right-hand panel), respectively, 
one before and one after tpeak- The mass distribution lies 
in the range 10^“^ M©, with few clumps as massive as the 
MBHs. These very massive clumps typically form after tpeak, 
most probably due to gas accretion from low-density regions 
and to mergers between less massive clumps, and eventu¬ 
ally will merge with the gas overdensity surrounding each 
MBHs. When one of these more massive clumps manages to 
approach an MBH at close range, then a transient MBH- 
clump binary system forms, strong gravitational perturba¬ 
tions develop, and the MBH orbit greatly deviates from 
its original path. This is the reason behind the ‘wiggling 
orbits’ seen in Fig. right-hand panel. The typical BH- 
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Figure 3. Face-on gas density map for run ‘ThFBl’ at time 
t = 2.1 Myr (i.e., when the star formation rate is maximum, 
see Fig. |^. The gas shocked after the first disc collision frag¬ 
ments into a large number of small clumps which very rapidly 
convert gas into new stellar mass. The black dots correspond to 
the positions of the two MBHs. 


clump distance when the transient binary system forms is 
~ 10 — 20 pc, which is always resolved with a number of 
cells > 10, thanks to the refinement prescription described 
in Lupi, Haardt & Dotti (2015), allowing us to accurately 


resolve the BH-clump close interaction. 

In the case of run ‘ThFBh’, because of the relatively 
higher density threshold for SF, a slightly larger number of 
more massive clumps forms, resulting in the more disturbed 
orbits (and faster decay) seen in Fig. 

We have then compared the above analysis regarding 
the MBH dynamics with runs employing the aforementioned 
blast wave-like feedback from SNe (BWFB-like runs). As 
discussed before, this feedback implementation aims at de¬ 
scribing non-thermal processes in the aftermath of SNa ex¬ 
plosions. We hnd that the dynamical evolution of the MBHs 
is largely independent upon the details of the SNa feed¬ 
back employed, making MBH dynamics results fairly robust 
against the different implementations of sub-grid physics. 


3.2 Gas dynamics 

We discuss here the dynamics of the gas during the merger 
event. We focus on the case with the low-density threshold 
for SF (run ‘ThFBF), keeping in mind that the higher den¬ 
sity case produces a qualitatively and quantitatively similar 
outcome. 

Fig.j^shows the gas distribution around the MBHB af¬ 
ter t = 11 Myr. On large scale (left-hand panel), the relic 
disc resulting from the collision of the progenitor discs is al¬ 
most totally disrupted because of SNa feedback. This resid¬ 
ual structure is counter-rotating relative to the MBHB or¬ 
bit. On scales of the order of few pc (right-hand panel), the 
gas which has not been converted into stellar particles set¬ 
tles in a circumbinary disc, with a total mass of few 10^ Mq. 
The small disc corotates with the MBHB thanks to the drag¬ 
ging of gas by the MBHs during their inspiral towards the 
centre. Note that this implies that the angular momentum 


Figure 4. Total mass in clumps for run ‘ThFBl’. The red dia¬ 
monds correspond to the times at which we computed the clump 
mass distribution shown in Fig. 


of the residual gas changed sign during the evolution of the 
system. 

We report in Fig. the evolution of the modulus of 
MBH orbital angular momentum and compared it to the 
modulus of the total angular momentum of the gas which 
is the closest to the MBHs in the simulation, dehned as 
the gas within a sphere of radius equal to 0.5 times the 
MBH separation. We observe that at the beginning of the 
simulation the angular momentum of the gas is larger than 
that of the MBHs, and we remind that the gas is counter¬ 
rotating. After > 4 Myr, the angular momentum associated 
with the MBH orbit exceeds that of the gas and in principle 
there are the conditions for a change in the sign of the gas 
angular momentum, being dragged by the MBHs. The gas 
angular momentum actually changes sign after ~ 9 Myr, 
when the MBH separation is ~ 45 pc. At this evolutionary 
stage, a large fraction (> 90%) of the initial gas mass is 
already converted in stellar particles. After ~ 10 Myr, when 
SNe start to explode, the released energy is radiated away 
by the small amount of residual gas, which is however unable 
to form further stellar mass at a comparable rate. In other 
words, star formation is not halted by SNa feedback, rather 
by gas consumption. 

Concerning the impact of blast wave feedback (BWFB- 
type runs), as expected it does not alter the gas dynamics 
for a time ~ AtsN (at that point the two MBHs have already 
reached the centre of the system). After that time, the al¬ 
most simultaneous SNa events release a fairly large amount 
of energy which heats the gas up but is not radiated away. 
The net result is that the remaining gas is pushed at very 
large distances from the MBHB (up to ~ 500 pc) by the 
increased pressure. The MBHB lives then in a very low- 
density environment, and no circumbinary disc is formed on 
any scale. 


3.3 Prompt SNa explosions 

Both the MBH and gas dynamics are unaffected by feedback 
for the hrst 10 Myr as this is the assumed lifetime of massive 
stars (and hence for the onset of SNa feedback). To test 
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Figure 1. MBH dynamical evolution for runs ‘ThFEP (solid black lines) and ‘ThFBh’ (dashed red lines). Left-hand panel: projected 
orbital evolution. Right-hand panel: MBH separation versus time. The blue dots correspond to the time of binary formation. 



Figure 2. Star formation in run ‘ThFBl’. Left-hand panel: total stellar mass (solid red line) in units of the initial disc mass and 
the residual gas mass in units of (dashed blue line) as a function of time. Right-hand panel: star formation rate versus time. 


how our results depend upon such choice, we cosnider the 
extreme case of AtsN = 0 Myr, i.e, massive stars explode as 
soon as they form. 

We find that, as long as the SNa feedback is governed 
by thermal processes, only small differences in the MBH 
dynamics exist compared to the standard delay case pre¬ 
viously discussed. This similarity occurs because the SNa 
energy is mostly released in high-density clumps, where gas 
cools down very rapidly, and the clumps can survive the ex¬ 
plosion. As a consequence, star formation can proceed until 
almost all clump gas is consumed. 

Large differences occur instead if, along the AtsN = 0 
assumption, we employ the blast wave recipe for SNa feed¬ 
back. In Fig. we compare the projected MBH orbits (left- 
hand panel) and the MBH separation versus time (right- 
hand panel) for runs BWFBh_prompt and ThFBh_prompt. 
In the case of blast wave like feedback, the orbital decay is 
slower, with a typical binary formation time-scale of > 13 
Myr. The difference is due to the early SNa explosions that, 
coupled with the blast wave-like feedback, tend to disrupt 


the gas clumps and to deplete the gas reservoir progressively 
forming around the MBHs. As a consequence, the two MBHs 
evolve in a lower density, smoother environment, where low- 
mass clumps are typically unable to induce strong orbital 
perturbations. The net result is a less disturbed orbital de- 
cay (Fig. l£l left-hand panel). 

We therefore conclude that in the case of prompt SNa 
explosions, contrary to the standard delay case, the dynam¬ 
ical evolution of the MBHs is strongly affected by the feed¬ 
back mechanism employed. The SF density threshold in¬ 
stead does not result in relevant differences anyway. 

While MBH dynamics is basically unaffected by the 
value of AtsN in the case of thermal SNa feedback, substan¬ 
tial differences occur in the dynamics of the gas component. 
Along with a small-scale corotating circumbinary disc, we 
do observe a further, much larger disc/ring-like structure 
on ~ 100 pc scale (see Figp^. Indeed, feedback from SNe 
does not occur suddenly after 10 Myr but it is instead di¬ 
luted in time, so that the (rapidly cooling) gas has time to 
readjust in a disc-like structure. Though several other possi- 
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Figure 5. Mass distribution of clumps in run ‘ThFEP. Left-hand panel: mass distribution at two selected times when the number of 
clumps is relatively small. Right-hand panel: same as left-hand panel, but at two times when the number of clumps is larger. The four 
selected times are marked as red dots in Fig.[^ 
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Figure 6. Face-on gas density maps for run ‘ThFBl’ around the MBHB at the end of the simulation (t ~ 11 Myr). Left-hand panel: 
on large scales the disc is almost totally disrupted because of SNa explosions. Right-hand panel: zoom in of the nuclear region where an 
inner corotating gas disc forms. 


ble explanations exist (e.g., secular evolution of the Galactic 
disc), it is tempting to associate such structure to the central 
molecular zone of the Milky Way ( Jones et al.||20lT ). It is 
interesting to note that the larger scale disc keeps memory 
of the initial angular momentum, and it is then counter¬ 
rotating with respect to the small inner circumbinary disc 
which is, as discussed above, dragged by the MBHB. 


The case of blast wave-like feedback is still different. We 
do not observe a disc-like structure, rather we find a mas¬ 
sive triaxial gas distribution surrounding the MBHB with 
density of few 10^ cm“^ (see Fig. 11). This difference is pro¬ 
duced by the different nature of the SNa feedback, which 
is in this case able to heat the gas and provide a pressure 
support large enough to prevent gas contraction. 


Because of the large fraction of gas available (due to 
the SNa feedback which reduces the net star formation by 
destroying gas clumps, as discussed above) the gas will con¬ 
tinue to cool down, resulting in alternated phases of star for¬ 
mation (due to gas cooling and contraction) and re-heating 
(due to SNa feedback). We observe a large number of dense 
gas streams flowing from low-density regions towards the 
centre where the MBHB resides. This large inflow will re¬ 
sult in a burst of star formation in the nucleus and in a 
following phase of SN explosions. The energy provided by 
SNe will then reheat the gas, stopping the contraction and 
eventually expand the entire gas structure into a less dense 
state. These alternated phases, if occurring for enough time, 
could convert a large fraction of gas into new stellar mass. 
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Figure 8. Same as 


Fig.but for runs ‘BWFBP (solid black lines) and ‘BWFBh’ (dashed red lines). 
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Figure 9. Same as Fig. but for runs ‘ThFBh.prompt’ (solid black lines) and ‘BWFBh.prompt’ (dashed red lines). The blue dots 
correspond to the time of binary formation. 


which could eventually form a massive nuclear stellar cluster 
surrounding the MBHB. 


4 DISCUSSION AND CONCLUSIONS 

By means of high-resolution, AMR hydrodynamical simu¬ 
lations, we explored the evolution of two massive gas discs 
hosting at their centre an MBH. The two discs are on an el¬ 
liptic orbit and merge, to mimic the encounter between two 
very gas-rich disc galaxies. To maximise the strength of the 
interaction, the orbital angular momentum is chosen to be 
antiparallel to the disc’s angular momenta. 

Strong shocks that develop during the merger of the 
two discs become sites of intense star formation, and stellar 
feedback alters significantly the thermal and dynamical state 
of the gas which undergoes a major transformation. Most 
of the gas is turned into new stellar particles through the 


formation of clumps of mass < lO^M©. Only few clumps 
form as massive as the MBHs, weighing lO^M©. 

We explored different SNa feedback recipes: the ther¬ 
mal and blast wave feedback, assuming a lifetime of ~ 10 
Myr for the massive stars. Furthermore, we considered a 
case in which prompt SNa explosion is coupled with both 
thermal and blast wave feedback. We find that the orbits of 
the two MBHs are perturbed due to their interaction with 
single clumps during the paring phase I, resulting in impul¬ 
sive kicks that imprint sudden changes in the direction and 
velocity of the orbit. Sinking times of ~ 10 — 20 Myr are 
found, considering the set of parameters used. The pairing 
phase terminates with the formation of a Keplerian binary. 


The MBH orbit is stochastic due to the presence of 
the gas clumps. However, we do not see a sizeable delay 
or spreading in the sinking time due to gas dumpiness, con¬ 
trary to what found in Fiacconi et al. (2013), where the level 
of stochasticity of the orbit was higher. We interpret this dif¬ 
ference as due to the geometry of the collision that mainly 
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Figure 10. Same as Fig. |^but for run ‘ThFBh.prompt’ at time t ~ 10 Myr. Left-hand panel: gas settles in a disc/ring like structure 
which is counter-rotating relative to the MBHs. Right-hand panel: zoom in of the region where an inner corotating gas disc forms around 
the MBHB visible on the east side of the left-hand panel. 
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Figure 11. Gas density maps for run ‘BWFBh.prompt’ around the MBHB at the end of the simulation (t 20 Myr). The gas settles in 
a triaxial structure with a denser central core. The core mass is ~ 10^ Mq within a radius ~ 25 pc. The upper panel shows the face-on 
view, while the two bottom panels show the edge-on views of the triaxial gas configuration. 


confines star formation along the oblique shock forming at 
the time of impact of the two discs, and to the fact that in 
our case the mass distribution of the clumps evolves as gas 
is turned into stellar mass which spread due to dynamical 
relaxation. Our simulated MBHs do not leave the orbital 
plane due to clump-induced kick, contrary to what seen in 


Roskar et al. (2015), as our simulation is strictly coplanar. 


We expect that an inclined encounter would lead to a change 
in the orbital plane also in our case, and this will be explored 
in future. During the pairing phase, the MBH dynamics is 


mostly affected by the presence of clumps and not by the 
recipe used to model the feedback processes. 

We note that, on the contrary, the gas distribution 
around the MBHs is significantly affected by feedback. Ther¬ 
mal feedback leaves no large-scale disc around the MBHB. 
Yet a residual corotating circumbinary disc of mass much 
smaller than the MBH mass forms around the two black 
holes which we expect will control the further spiral-in via 
migration-like mechanisms. 

Blast wave feedback is a way to model the expansion 
of SNa-driven bubbles. With the code it is then possible to 
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Figure 7. Evolution of the moduli of the angular momenta for 
run ‘ThFBl’. The MBHB total angular momentum (solid black 
line) is compared with the gas angular momentum (dashed red 
line), defined as the total angular momentum of the gas within a 
sphere centred on one MBH whose radius equals half the MBHB 
separation. Both curves have been normalised to the initial an¬ 
gular momentum of the MBHB. 


mimic the ballistic phase of the shock triggered by the SNa 
explosion. As cooling is shut off in this phase, a multiphase 
gas forms and the sweeping of the gas induced by the blast 
wave leads to the almost complete evacuation of gas. The 
MBHB thus inhabits a region completely devoid of gas. Blast 
wave feedback in the prompt scenario leads instead to a 
configuration in which the MBHB is surrounded by a gas 
cloud with little angular momentum and triaxial in shape. 

The lesson to learn is that star formation in merging 
gaseous discs is a key process which affects the physical state 
of the gas in the surroundings of the MBHs. Under these cir¬ 
cumstances it is difficult to predict the actual distribution 
of gas when the most active phase of the merger has sub¬ 
sided, as the outcome depends upon the modelling and on 
sub-grid physics, and firm conclusions should be taken with 
caution. Still, the presence of cool gas has deep implications 
for the evolution and observability of close MBHBs. First, 
the evolution of a binary on sub-pc scales towards the co¬ 
alescence is strongly dependent on the gaseous and stellar 


distribution in its immediate surroundings (Colpi & Dotti 


2011 for a review). The time-scale of the MBHs shrinking 
on sub-pc scales is of fundamental importance as it affects 
the expected rate of binaries possibly observable as gravi¬ 
tational wave sources. This is particularly true in mergers 
between gas-rich galaxies, a fraction of which can host bi¬ 
naries detectable by future space-based gravitational wave 
detectors such as eLISA ( Amaro-Seoane et al.||2013 ). Sec¬ 
ondly, the presence of gas is a necessary condition for the 
possible detection of the binary during the hardening phase 
(Dotti, Sesana & Decarli 2012) as well as for pinpointing 
an electromagnetic counterpart of the MBHB coalescence 
(see e.g. Schnittman 2013 Bogdanovic 2015). The lack of 
a clear consensus on the processes shaping the environment 
of MBHBs, whose evolution actually depends on the physi¬ 
cal modelling, and the lack of observations available on the 
small scales we considered in our work witness the need of 
investigating a wider range of parameters. 


We thank R. Teyssier and L. Paredi for many fruitful discus¬ 
sions. We acknowledge financial support from Italian MIUR, 
through PRIN 2010-2011. Simulations were run on the EU- 
RORA cluster at CINECA and on the Lucia cluster at 
DiSAT, University of Insubria. 
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